function y = transitions(xvals,Z,Tstep)  

%This function uses hard coded values for the locations of preference clusters and
%their weights, and computes the elements of the switching probability
%matrix.  There are two options, either switching probabilities are
%conditionally uniformly distributed, or they depend on the Mahalanobis
%distance as described in the appendix of the paper.  The output is a
%matrix whose first column is the pure discount factor (scaled to 10 year time step), second column is
%values of eta, and remaining columns are the elements of the transition
%matrix.

flag = 0;   %Change flag to 1 to use Mahalananobis distance when computing transition matrix.

load('deltaeta_data.mat');

Ddata = deltaeta_data;

%First column of C below is the 10-year beta, second column is eta.

C = [0.48519	0.85
0.98341	1.85
0.92036	4.6667
0.55839	4
0.86211	1.8421
0.65012	0.74
0.94453	2.9286
0.96573	1.037
0.79849	0.83562
0.95739	0.4053];

%Cluster weights

weight = [0.011565475764550
   0.219650388150392
   0.017333016188854
   0.011565475764552
   0.109825575046945
   0.028895265388023
   0.040456505023887
   0.265894959777211
   0.167632921791943
   0.127178937797554]';



%This little routine computes a distance matrix for the clusters, based on Mahalanobis
%distance
S = cov(Ddata(:,1:2));
invS = inv(S);
D = zeros(Z);
for i = 1:Z
    for j = 1:Z
        D(i,j) = sqrt((C(i,:)-C(j,:))*invS*(C(i,:)-C(j,:))');  %Mahalanobis distance.
    end
end
Dmax = max(max(D));
lambda = ((0.1)^(-1)-1)/Dmax;  %The parameter lambda is chosen so that the biggest distance between preferences has a weighting of 0.1 when we compute the switching probabilities.

V = [];


%Compute the matrix of switching probabilities (annual) for each value of
%x (non-dogmatism parameter)
for x = xvals
    
    V = [V; x*ones(1,Z+2)];

    M = zeros(Z);

    for i = 1:Z
        A = (sum(weight./(1+lambda*D(i,:))))^(-1);  %This is a normalization constant for the switching probabilities for type i, used in the Mahalnonobis case.  It ensures that sum of all switching probabilities for type i is 1.
        for j = 1:Z
            if j==i
                
                if flag == 0
                    
                    M(i,j) = (1-x + x*weight(j));
                    
                elseif flag == 1
                
                    M(i,j) = 1-x + x*A*weight(j);
                end
                
                %
            else
                if flag == 0
                    M(i,j) = x*weight(j);
                elseif flag == 1 
                    M(i,j) = x*A*weight(j)*(1/(1+lambda*D(i,j)));
                end
            end
        end    
    end
    
    M = M^(Tstep);    %Scale transition matrix to time step.
    
    V = [V; C(:,1), C(:,2), M];

end

y = V;
%



